%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%{
Article citation:
F. Jord�n-Cuebas, U. Krogmann, Laundry Energy Consumption in Multistory Buildings: Technology versus Laundering Practices, Journal of Building Engineering, https://doi.org/10.1016/j.jobe.2026.115189.

Source code/data citation:
F. Jord�n-Cuebas, U. Krogmann, Replication Data for: Laundry Energy Consumption in Multistory Buildings: Technology versus Laundering Practices, Harvard Dataverse, https://doi.org/10.7910/DVN/BTHCFL.
%}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

clc
clearvars
load('MCNA_inputdata_v1')
g = 'one';
b = input.(g);

for mcna_m = [1 0]
	for mcna_s = [0 1 2]
	
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	%%% USER INPUT - SELECT MODEL AND SENSITIVITY ANALYSIS
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	whichmodel = mcna_m; % SELECT MODEL: 1 = based on literature input data; 0 = based on case study building + literature input data
	whichsensi = mcna_s; % SELECT SENSITIVITY ANALYSIS (SI = sensitivity index): 0 = none; 1 = SImax; 2 = SImin
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	
	% Sensitivity index parameters
	if whichsensi == 1
		senselect = input.(g).mu(184:187)'; % SImax
	elseif  whichsensi == 2
    	senselect = input.(g).min(184:187)'; % SImin
	else
    	senselect = input.(g).mass(184:187)';% none
	end
	sensicalc = senselect(1);
	senseindex = senselect(2);
	sensitype = senselect(3);
	sensichange = senselect(4);
	redist = 0;
	sensibase = 1;
	
	for scenario = [s_nameSC.Var1]' % Scenario names
		
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		%%% Input data
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		n = input.(g).name(1:4)'; % Fabric program names
		lb = input.(g).name(123:130)'; % Load mass names
		x = input.(g).name(5:8)'; % Wash water temperature names
		y = input.(g).name(9:11)'; % Dryness level names
		e = input.(g).name(33:33)'; % Water heater names
		hx = input.(g).name(86:90)'; % Washing frequency of use names
		hxn = input.(g).mass(86:90)'; % Load mass household size correction factor
		n2 = input.(g).name(236:243)'; % RMC by program and RPM names
		n3 = input.(g).name(269:276)'; % RMC names
		n3mass = input.(g).mass(269:277)'; % Dryer electricity consumption, load mass
		n3rmc = input.(g).mu(269:276)'; % RMCnet, correction factor for dryer electricity consumption
		n3elect = input.(g).min(269:276)'; % Dryer electricity consumption, electricity consumption
		nhot = input.(g).name(310:312)'; % Fraction of main cycle water consumption that is building hot water for hot water names 
		nwarm = input.(g).name(313:315)'; % Fraction of main cycle water consumption that is building hot water for warm water names 
		n2r = input.(g).mass(236:243)'; % RPM bins
		ss = 10000; % Number of iterations in Monte Carlo simulation
		
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		%%% Triangular distributions
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		rng('default');
		triangle = @(table, row) random(makedist('Triangular','a',b{row,4},'b',b{row,3},'c',b{row,5}),ss,1);
		
		for i = 1:length(n)
			mass.(n{i}).s = triangle(b, 36+(i-1)); % Load mass by fabric program; replace line with mass.(n{i}).s = triangle(b, 212+(i-1)) to simulate widening the load mass data range to that of washing frequency
			massfull.(n{i}).s = random(makedist('Uniform','Lower',input.(g).min(67),'Upper',input.(g).max(67)),ss,1).*input.(g).mu(67); % Scenario #15 - Large loads
			masssmall.(n{i}).s = 1.1; % Scenario #2 - Small loads
			rinse.extra.(n{i}).s = triangle(b, 44); % Fraction of loads with an extra rinse cycle
			rinse.extraall.(n{i}).s = triangle(b, 210); % Scenario #11 - Extra rinse
		end
		
		for i = 1:length(n)
			for kg = 1:length(lb)
        		waterbin.(n{i}).(lb{kg}) = triangle(b, (123+((i-1)*8))+(kg-1)); % Washer water consumption in the main and rinse cycles
        		electbin.(n{i}).(lb{kg}) = triangle(b, (156+(kg-1))); % Washer operational electricity consumption
				labwaterbin.(n{i}).(lb{kg}) = triangle(b, (278+((i-1)*8))+(kg-1)); % Washer water consumption in the main and rinse cycles in the lab washers
			end
		end
		
		for i = 1:length(n)
			for i2 = 1:length(n2)
        		moistbin.(n{i}).(n2{i2}) = triangle(b, (236+((i-1)*8))+(i2-1)); % RMC of the wet load by program
			end
		end
		
		for i = 1:length(n)
    		topload.(n{i}) = triangle(b,189); % Scenario #1 - Conventional top-loading washer water consumption
			efftopload.(n{i}) = triangle(b,190); % Scenario #3 - Efficient top-loading washer water consumption
		end
		
		zr1=1;
		
		for s = 1:length(x)
			tubtemp2.(x{s}).s = triangle(b, 45+(s-1)); % Washer wash water temperature in the main cycle
		end
		
		for s = 1:length(x)
    		if whichmodel ==1
				if scenario == 4
					tempmain.(x{s}) = triangle(b, 112+(s-1)); % Scenario #17 - Cold wash water temperature
				elseif scenario == 5
					tempmain.(x{s}) = triangle(b, 173+(s-1)); % Scenario #5 - Highest appropriate washer temperature setting
				else    
					tempmain.(x{s}) = triangle(b, 101+(s-1)); % Fraction of loads by main cycle water temperature
				end 
    		else
        		tempmain.(x{s}) = triangle(b, 5+(s-1)); % Fraction of loads by main cycle water temperature for the case study building
    		end    
		end
		
		for s = 1:length(x)
			tubtempsupp.(x{s}).s = triangle(b, 91+(s-1)); % Washer supplemental wash water temperature in the main cycle
		end
		
		for i = 1:length(n)
			if whichmodel == 1
				if scenario == 3
					fracpgrm.(n{i})= triangle(b, 194+(i-1)); % Scenario #8 - Permanent press program only
				elseif scenario == 23
					fracpgrm.(n{i})= triangle(b, 206+(i-1)); % Scenario #9 - Cotton program only				
				elseif scenario == 25
					fracpgrm.(n{i})= triangle(b, 202+(i-1)); % Scenario #16 - Hand wash program only		
				else	
           		fracpgrm.(n{i})= triangle(b, 1+(i-1)); % Fraction of loads by program
				end
			else
				fracpgrm.(n{i})= triangle(b, 169+(i-1)); % Fraction of loads by program for the case study building
			end
		end
		
			if scenario == 30
				for i = 1:length(n)
					mainrpm.(n{i})= min(triangle(b, 29+(i-1))); % Scenario #4 - Low final spin speed
				end
			elseif scenario == 31 
				for i = 1:length(n)
					mainrpm.(n{i})= max(triangle(b, 29+(i-1))); % Scenario #18 - New high spin speed washer
				end
			elseif scenario == 13 
				for i = 1:length(n)
					mainrpm.(n{i})= triangle(b, 316); % Scenario #1 - Conventional top-loading washer
				end
			elseif scenario == 16
				for i = 1:length(n)
					mainrpm.(n{i})= triangle(b, 317); % Scenario #3 - Efficient top-loading washer
				end
			else
				for i = 1:length(n)
					mainrpm.(n{i})= triangle(b, 29+(i-1)); % Washer final spin speeds by program
				end
			end
		
		for u = 1:length(y)
			rema.dryer.lit.(y{u}) = triangle(b, 57+(u-1)); % RMC of the dry load by dryness level
		end
		
		for u = 1:length(y)
    		if scenario == 9
        		drymain.(y{u})(:,1) = triangle(b, 181+(u-1)); % Scenario #6 - Higher dryness level in dryer
    		elseif scenario == 10
        		drymain.(y{u})(:,1) = triangle(b, 120+(u-1)); % Scenario #19 - Lower dryness level in dryer
    		else    
        		drymain.(y{u}) = triangle(b, 9+(u-1)); % Fraction of loads by dryness level
    		end 
		end
		
		for w = 1:length(e)
			dhw.(e{w}) = 1./triangle(b, 33); % Building hot water system losses factor
		end
		
		bldghhprecat = triangle(b, 168); % Household size for case study building
		autonormalfactor = triangle(b, 105); % Dryer automatic termination correction factor for normal dry
		autoextrafactor = triangle(b, 106); % Dryer automatic termination correction factor for more dry
		
		recshh = triangle(b, 167); % Household size
		
		for htx = 1:length(hx)
			recs2018(~recs2018.wshloadcapd,:) = [];
			hh1 = recs2018.wshloadcapd; % Washing frequency based on RECS 2018
			hh2 = recs2018.NHSLDMEM2; % Household size
			hhs = htx;
			[distr] = extdist(hh1, hh2, hhs, ss);
			bldghhx.(hx{htx}) = distr; % Washing frequency by household size
		end
		
		storage_temp = triangle(b, 68); % Hot water temperature in the apartments in case study building
		bldg_temp = triangle(b, 69); % Water heaters temperature setpoint in the case study building
		
		if scenario == 22
			source.electricity.s = triangle(b, 223); % Scenario #13 - Increasing contribution of renewable energy sources
		else
			source.electricity.s = triangle(b, 222); % Source-to-site energy factor for electricity
		end
		
		source.gas.s = triangle(b, 65); % Source-to-site energy factor for natural gas
		mains_temp = tubtemp2.cold.s; % Water temperature in the city mains
		
		fsg = triangle(b, 66); % Fraction of main and rinse cycles water consumption in the main cycle
		fsgmpre = triangle(b, 95); % Fraction of main and rinse cycles water consumption in the extra rinse cycle
		
		if scenario == 15
    		dryfract = triangle(b, 70)./2; % Scenario #21 - Line drying half of the clothes   
		else
    		dryfract = triangle(b, 70); % Fraction of loads dried in the dryer
		end
		
		heatpumpenergysavings = 1 - triangle(b,188); % Scenario #20 - Heat pump clothes dryer
		
		spraytechelect = triangle(b, 191); % Scenario #14 - Front-loading washer with spray technology energy consumption
		spraytechwater = triangle(b, 192); % Scenario #14 - Front-loading washer with spray technology water consumption
		
		moreclothing = triangle(b, 211); % Scenario #7 - Increasing amount of clothes by 5 %; replace line with moreclothing = triangle(b, 211).*10 to simulate increasing amount of clothes 10x = 50 %
		
		for i = 1:length(n)
			rmcwetlab.(n{i}) = triangle(b, 228+(i-1)); % RMC of the wet load by fabric program for the case study building
		end
		
		for ihot = 1:length(nhot)
			hotwatfrac.(nhot{ihot}) = triangle(b, 310+(ihot-1)); % Faction of main cycle water consumption that is building hot for the hot cycle in lab washer
		end
		
		for iwarm = 1:length(nwarm)
			hotwatfrac.(nwarm{iwarm}) = triangle(b, 313+(iwarm-1)); % Faction of main cycle water consumption that is building hot for the warm cycle in lab washer
		end
		
		n2corr = 1 - triangle(b, 268); % Bone dry load mass correction
		
		func = @nanmedian; % Function handle for median
		func25 = @(data) quantile(data,0.25); % Function handle for 25th percentile
		func75 = @(data) quantile(data,0.75); % Function handle for 75th percentile
		
		outlier = 'kswo';
		
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		% Model, sensitivity analysis and scenarios
		%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
		sensi = 72; % Number of variables in the sensitivity analysis
		for snt = 0:sensi % Initial default is 0:sensi
			clc
			
			name = 'MCNA';   
			fprintf('%s at work...',name);
			
			for sn = 1:sensi
    			eval(['s_',num2str(sn),'=',num2str(sensibase),';']);
			end	
				
			if snt == 0
    			for del = 1:sensi
        			eval(['s_',num2str(del),'=',num2str(sensibase),';']);
    			end	
			else
    			if 	eval(['s_',num2str(sn)]) ~= eval(['s_',num2str(snt)]);
        			eval(['s_',num2str(snt),'=',num2str(sensibase),';']);
    			else
        			eval(['s_',num2str(snt),'=',num2str(sensichange),';']);
    			end
			end
			
			%% Fraction of loads by program sensitivity analysis
			temp_program = [];
			for i = 1:length(n) 
    			temp_program = sum(horzcat(fracpgrm.(n{i}), temp_program),2); % Fraction of loads by program input from triangular distribution
			end
			
			for i = 1:length(n)
    			tempprogramstd.(n{i}) = fracpgrm.(n{i})./temp_program; % Standardization
			end
			
			for i = 1:length(n)
    			sntset = 9:12; % Variables for sensitivity analysis
    			absval = std(tempprogramstd.(n{i})); % Standard deviation
    			[ndiff, setdi, ingroup] = preppgrm (i, n, sntset); % Function preparation
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist); % Function for sensitivity analysis
    			if senseindex == 0
        			frac.(n{i})(:,1) = c2*tempprogramstd.(n{i}) + d2;
    			else
					if c2 > 1
						frac.(n{i})(:,1) = max(tempprogramstd.(n{i})); % Sensitivity index (SI) maximum value 
					elseif c2 < 1
						frac.(n{i})(:,1) = min(tempprogramstd.(n{i})); % Sensitivity index (SI) minimum value
					else c2 = 1;
						frac.(n{i})(:,1) = tempprogramstd.(n{i});
					end
    			end
			end
			
			temp_frac = []; 
			for i = 1:length(n)
				temp_frac = sum(horzcat(frac.(n{i}), temp_frac),2); % Fraction of loads by program modeled
			end
			
			% NOTE: Comments will be added only to new lines from herein.
			
			%% Fraction of loads by main cycle water temperature sensitivity analysis
			temp_main = [];
			for s = 1:length(x)
    			temp_main = sum(horzcat(tempmain.(x{s}), temp_main),2);
			end
			
			for s = 1:length(x)
    			tempmainstd.(x{s}) = tempmain.(x{s})./temp_main;
			end
			
			for s = 1:length(x)
    			sntset = 14:17;
    			absval = std(tempmainstd.(x{s}));
    			[ndiff, setdi, ingroup] = preptemp (s, x, sntset);
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
    			if senseindex == 0
        			washtempfract.(x{s})(:,1) = c2*tempmainstd.(x{s}) + d2;
				else
					if c2 > 1
						washtempfract.(x{s})(:,1) = max(tempmainstd.(x{s}));
					elseif c2 < 1
						washtempfract.(x{s})(:,1) = min(tempmainstd.(x{s}));
					else c2 = 1;
						washtempfract.(x{s})(:,1) = tempmainstd.(x{s});
					end
    			end
			end
			
			%% Household size sensitivity analysis
			if whichmodel == 1 % If condition (whichmodel == 1) met then literature input data model
				bldghh(:,1) = recshh;
				for m5 = 1:length(zr1)          
					sntset = 71:71;
					absval = std(recshh);
					[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
					[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
					if senseindex == 0
						bldghh(:,1) = c2*recshh + d2; 
					else
						if c2 > 1
							bldghh(:,1) = max(recshh);
						elseif c2 < 1
							bldghh(:,1) = min(recshh);
						else c2 = 1;
							bldghh(:,1) = recshh;
						end
					end
				end
			
			else % If condition not met then case study building + literature input data model
				bldghh(:,1) = bldghhprecat;
			end
			
			%% Load mass corrected by household size
			for htx = 1:length(hx)
				bldghhxbinmass.(hx{htx}) = ((bldghh >= htx) & (bldghh < (htx)+1)).*(hxn(htx));
			end
			
			temp_hhsizeloadmass = [];
			for htx = 1:length(hx) 
				temp_hhsizeloadmass = sum(horzcat(bldghhxbinmass.(hx{htx}), temp_hhsizeloadmass),2);
			end
			
			%% Washing frequency by household size sensitivity analysis
			for htx = 1:length(hx)
    			sntset = 20:24;
    			absval = std(bldghhx.(hx{htx}));
    			[ndiff, setdi, ingroup] = prephhs (htx, hx, sntset);   
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);   
    			if senseindex == 0
        			bldghhx2.(hx{htx})(:,1) = c2*bldghhx.(hx{htx}) + d2;
				else
					if c2 > 1
						bldghhx2.(hx{htx})(:,1) = max(bldghhx.(hx{htx}));
					elseif c2 < 1
						bldghhx2.(hx{htx})(:,1) = min(bldghhx.(hx{htx}));
					else c2 = 1;
						bldghhx2.(hx{htx})(:,1) = bldghhx.(hx{htx});
					end
    			end 
			end
			
			for htx = 1:length(hx)
    			bldghhxbinload.(hx{htx}) = ((bldghh >= htx) & (bldghh < (htx)+1)).*bldghhx2.(hx{htx});
			end
			
			temp_hhsizeload = [];
			for htx = 1:length(hxn)
    			temp_hhsizeload = sum(horzcat(bldghhxbinload.(hx{htx}), temp_hhsizeload),2);
			end
			
			for i = 1:length(n) 
        			freq.bldg2.(n{i}) = temp_hhsizeload.*frac.(n{i}); % Washing frequency by fabric program
			end
			
			%% Clothing
			for i = 1:length(n)
				clothing.(n{i}) = freq.bldg2.(n{i}).*(mass.(n{i}).s.*temp_hhsizeloadmass);
			end
			
			for i = 1:length(n)
				if scenario == 17 % Scenario #7 - Increasing amount of clothes by 5 %
					clothing.(n{i}) = freq.bldg2.(n{i}).*(mass.(n{i}).s.*temp_hhsizeloadmass).*(1+moreclothing); % Clothing
					freq.bldgb.(n{i}) = clothing.(n{i})./(mass.(n{i}).s.*temp_hhsizeloadmass); % Washing frequency
					mass.(n{i}).b = clothing.(n{i})./freq.bldgb.(n{i}); % Load mass
				elseif scenario == 23 | scenario == 3 | scenario == 25 % Scenario #9 - Cotton program only; Scenario #8 - Permanent press program only; and Scenario #16 - Hand wash program only
					clothing.(n{i}) = clothing_baseline.*frac.(n{i});
					freq.bldgb.(n{i}) = (clothing_baseline./(mass.(n{i}).s.*temp_hhsizeloadmass)).*frac.(n{i});
					mass.(n{i}).b = (clothing.(n{i})./freq.bldgb.(n{i})).*frac.(n{i});
				elseif scenario == 6 % Scenario #15 - Large loads
					mass.(n{i}).b =  massfull.(n{i}).s.*temp_hhsizeloadmass;
					freq.bldgb.(n{i}) = clothing.(n{i})./mass.(n{i}).b; 
				elseif scenario == 7 % Scenario #2 - Small loads
					mass.(n{i}).b = masssmall.(n{i}).s;
					freq.bldgb.(n{i}) = clothing.(n{i})./mass.(n{i}).b;
				else
					mass.(n{i}).b = clothing.(n{i})./freq.bldg2.(n{i});
					freq.bldgb.(n{i}) = clothing.(n{i})./mass.(n{i}).b;
				end
			end
			
			clothing_all = [];
			for i = 1:length(n)
				clothing_all = sum(horzcat(clothing.(n{i}), clothing_all),2);
			end
			
			for i = 1:length(n)
				freq.cap.(n{i})(:,1) = c2*freq.bldgb.(n{i});
			end
			
			%% Load mass by fabric program sensitivity analysis
			for i = 1:length(n)
    			sntset = 1:4;
				mass.(n{i}).bd = mass.(n{i}).b.*n2corr;
    			absval = std(mass.(n{i}).bd);
    			[ndiff, setdi, ingroup] = preppgrm (i, n, sntset);
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
    			if senseindex == 0
        			mass.washer.(n{i})(:,1) = c2*mass.(n{i}).bd + d2;
    			else
					if c2 > 1
						mass.washer.(n{i})(:,1) = max(mass.(n{i}).bd);
					elseif c2 < 1
						mass.washer.(n{i})(:,1) = min(mass.(n{i}).bd);
					else c2 = 1;
						mass.washer.(n{i})(:,1) = mass.(n{i}).bd;
					end
    			end
			end
			
			temp_mass = [];
			for i = 1:length(n)
				temp_mass = sum(horzcat(mass.washer.(n{i}).*(1./n2corr).*frac.(n{i}), temp_mass),2);
			end
			
			%% Washer water consumption in the main and rinse cycles by program and load mass sensitivity analysis
			for i = 1:length(n)
				for kg = 1:length(lb)
        			loadbinpre.(n{i}).(lb{kg}) = ((mass.washer.(n{i}) >= kg) & (mass.washer.(n{i}) < (kg+1))).*kg.*waterbin.(n{i}).(lb{kg});
				end
				
				temp_loadbinpre.(n{i}) = [];
				for kg = 1:length(lb)
        			temp_loadbinpre.(n{i}) = sum(horzcat(loadbinpre.(n{i}).(lb{kg}), temp_loadbinpre.(n{i})),2);
				end
				
				sntset = 35:38;
    			absval = std(temp_loadbinpre.(n{i}));
    			[ndiff, setdi, ingroup] = preppgrm (i, n, sntset);
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
				
    			if senseindex == 0
        			water.eff.(n{i})(:,1) = c2*temp_loadbinpre.(n{i}) + d2;
    			else
					if c2 > 1
						water.eff.(n{i})(:,1) = max(temp_loadbinpre.(n{i}));
					elseif c2 < 1
						water.eff.(n{i})(:,1) = min(temp_loadbinpre.(n{i}));
					else c2 = 1;
						water.eff.(n{i})(:,1) = temp_loadbinpre.(n{i});
					end
    			end
			end    
			
			for i = 1:length(n)
				for kg = 1:length(lb)
        			labloadbinpre.(n{i}).(lb{kg}) = ((mass.washer.(n{i}) >= kg) & (mass.washer.(n{i}) < (kg+1))).*kg.*labwaterbin.(n{i}).(lb{kg});
				end
					water.lab.tot.(n{i}) = [];
				for kg = 1:length(lb)
        			water.lab.tot.(n{i}) = sum(horzcat(labloadbinpre.(n{i}).(lb{kg}), water.lab.tot.(n{i})),2);
				end	
			end
			
			if whichmodel == 1
				if scenario == 13
					for i = 1:length(n)
						water.washer.(n{i}) = topload.(n{i}).*mass.washer.(n{i}); % Scenario #1 - Conventional top-loading washer
					end
				elseif scenario == 16
					for i = 1:length(n)
						water.washer.(n{i}) = efftopload.(n{i}).*mass.washer.(n{i});  % Scenario #3 - Efficient top-loading washer
					end
				elseif scenario == 14
					for i = 1:length(n)
						water.washer.(n{i}) = water.eff.(n{i}).*(1-spraytechwater); % Scenario #14 - Front-loading washer with spray technology
					end			
				else
					for i = 1:length(n)		
						water.washer.(n{i}) = water.eff.(n{i}); % Per load washer water consumption	
					end
				end
			
				temp_water = [];
				for i = 1:length(n)
					temp_water = sum(horzcat(water.washer.(n{i}).*freq.cap.(n{i}), temp_water),2); % Daily per capita washer water consumption	
				end
			else
				for i = 1:length(n)
        			water.washer.(n{i}) = water.lab.tot.(n{i});  % Per load washer water consumption for case study bulding
				end
				
				temp_water = [];
    			for i = 1:length(n)
        			temp_water = sum(horzcat(water.washer.(n{i}).*freq.cap.(n{i}), temp_water),2); % Per capita daily washer water consumption for case study building
    			end	
			end
			
			temp_water_kg = [];
			for i = 1:length(n)
				temp_water_kg = sum(horzcat((water.washer.(n{i})./mass.washer.(n{i})(:,1)).*frac.(n{i}), temp_water_kg),2); % Specific washer water consumption
			end	
			
			temp_water_ld = [];
			for i = 1:length(n)
				temp_water_ld = sum(horzcat(water.washer.(n{i}).*frac.(n{i}), temp_water_ld),2); % Per load washer water consumption
			end	
			
			%% Washer hot water consumption
			for s = 1:length(x)
    			tubtemp.(x{s}).s = tubtemp2.(x{s}).s; 
			end
			
			for i = 1:length(n)
				for s = 1:length(x)
        			for m5 = 1:length(zr1) % Fraction of main cycle water consumption that is building hot water sensitivity analysis  
            			sntset = 19:19;
            			absval = std((tubtemp.(x{s}).s - mains_temp)./(storage_temp - mains_temp));
            			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
            			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
						if senseindex == 0
							tempihot.(x{s})(:,1) = c2*((tubtemp.(x{s}).s - mains_temp)./(storage_temp - mains_temp)) + d2;
						else
							if c2 > 1
								tempihot.(x{s})(:,1) = max((tubtemp.(x{s}).s - mains_temp)./(storage_temp - mains_temp));
							elseif c2 < 1
								tempihot.(x{s})(:,1) = min((tubtemp.(x{s}).s - mains_temp)./(storage_temp - mains_temp));
							else c2 = 1;
								tempihot.(x{s})(:,1) = (tubtemp.(x{s}).s - mains_temp)./(storage_temp - mains_temp);
							end
						end
        			end       
 			
        			for m5 = 1:length(zr1) % Fraction of main and rinse cycles water consumption in the main cycle sensitivity analysis
            			sntset = 18:18;
            			absval = std(fsg);
            			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
            			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
						if senseindex == 0
							hotwater.washer.(n{i}).(x{s})(:,1) = (c2*fsg + d2).*water.washer.(n{i}).*tempihot.(x{s}).*washtempfract.(x{s});
						else
							if c2 > 1
								hotwater.washer.(n{i}).(x{s})(:,1) = max(fsg).*water.washer.(n{i}).*tempihot.(x{s}).*washtempfract.(x{s});
							elseif c2 < 1
								hotwater.washer.(n{i}).(x{s})(:,1) = min(fsg).*water.washer.(n{i}).*tempihot.(x{s}).*washtempfract.(x{s});
							else c2 = 1;
								hotwater.washer.(n{i}).(x{s})(:,1) = fsg.*water.washer.(n{i}).*tempihot.(x{s}).*washtempfract.(x{s});
							end
						end
        			end
				end
			end
			
			for i = 1:length(n) 
    			temp_hotwash.(n{i}) = [];
				for s = 1:length(x)
        			temp_hotwash.(n{i}) = sum(horzcat(hotwater.washer.(n{i}).(x{s}), temp_hotwash.(n{i})),2); % Per load washer hot water consumption
				end
			end
			
			temp_hotwashf_lit = [];
			for i = 1:length(n)
    			temp_hotwashf_lit = sum(horzcat(temp_hotwash.(n{i}).*freq.cap.(n{i}),temp_hotwashf_lit),2); % Per capita daily washer hot water consumption
			end
			
			for i = 1:2 %  Washer hot water consumption for cotton and permanent press programs for the case study building
    			for s = 1:2 % Washer hot water consumption for xhot and hot wash water temperatures for the case study building
					if mass.washer.(n{i}) < 3
						water.lab.hot.(n{i}).(x{s}) = water.lab.tot.(n{i}).*hotwatfrac.(nhot{1}).*washtempfract.(x{s});
					elseif mass.washer.(n{i}) >= 4
						water.lab.hot.(n{i}).(x{s}) = water.lab.tot.(n{i}).*hotwatfrac.(nhot{3}).*washtempfract.(x{s});
					else
						water.lab.hot.(n{i}).(x{s}) = water.lab.tot.(n{i}).*hotwatfrac.(nhot{2}).*washtempfract.(x{s});
					end
    			end
			end
			
			for i = 1:length(n) % Washer hot water consumption for all programs for the case study building
    			for s = 3:3 % Washer hot water consumption for warm wash water temperature for the case study building
					if mass.washer.(n{i}) < 3
						water.lab.hot.(n{i}).(x{s}) = water.lab.tot.(n{i}).*hotwatfrac.(nwarm{1}).*washtempfract.(x{s});
					elseif mass.washer.(n{i}) >= 4
						water.lab.hot.(n{i}).(x{s}) = water.lab.tot.(n{i}).*hotwatfrac.(nwarm{3}).*washtempfract.(x{s});
					else
					water.lab.hot.(n{i}).(x{s}) = water.lab.tot.(n{i}).*hotwatfrac.(nwarm{2}).*washtempfract.(x{s});
					end		
    			end
			end
			
			for i = 1:2
    			temp_water_hot.(n{i}) = [];
    			for s = 1:3
        			temp_water_hot.(n{i}) = sum(horzcat(water.lab.hot.(n{i}).(x{s}), temp_water_hot.(n{i})),2); % Per load washer hot water consumption for cotton and permanent press programs and xhot and hot wash water temperatures for the case study building
    			end
			end
			
			for i = 1:length(n)
    			temp_water_hot.(n{i}) = [];
    			for s = 3:3
        			temp_water_hot.(n{i}) = sum(horzcat(water.lab.hot.(n{i}).(x{s}), temp_water_hot.(n{i})),2); % Per load washer hot water consumption for all programs and warm wash water temperature for the case study building
    			end
			end
			
			temp_water_hotlab = [];
			for i = 1:length(n)
    			temp_water_hotlab = sum(horzcat(temp_water_hot.(n{i}).*freq.cap.(n{i}), temp_water_hotlab),2); % Per capita daily washer hot water consumption for the case study building
			end
			
			if whichmodel == 1 
    			temp_hotwashf = temp_hotwashf_lit;
			else     
    			temp_hotwashf = temp_water_hotlab;
			end
	
			%% Washer extra rinse water consumption
			if whichmodel == 1
				for m5 = 1:length(zr1) % Fraction of main and rinse cycles water consumption in the extra rinse cycle sensitivity analysis
					sntset = 30:30;
					absval = std(fsgmpre);
					[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
					[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
					if senseindex == 0
					fsgm(:,1) = c2*fsgmpre + d2;
					else
						if c2 > 1
							fsgm(:,1) = max(fsgmpre);
						elseif c2 < 1
							fsgm(:,1) = min(fsgmpre);
						else c2 = 1;
							fsgm(:,1) = fsgmpre;
						end
					end
				end
	
				for i = 1:length(n)
					sntset = 54:57; % Fraction of loads with an extra rinse cycle sensitivity analysis
					absval = std(rinse.extra.(n{i}).s);
					[ndiff, setdi, ingroup] = preppgrm (i, n, sntset);
					[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
					if senseindex == 0
						extrarinse.washer.(n{i})(:,1) = c2.*rinse.extra.(n{i}).s + d2;
					else
						if c2 > 1
							extrarinse.washer.(n{i})(:,1) = max(rinse.extra.(n{i}).s);
						elseif c2 < 1
							extrarinse.washer.(n{i})(:,1) = min(rinse.extra.(n{i}).s);
						else c2 = 1;
							extrarinse.washer.(n{i})(:,1) = rinse.extra.(n{i}).s;
						end
					end
				end
				if scenario == 26 % Scenario #11 - Extra rinse
					temp_extrarinse = [];
					for i = 1:length(n)
						temp_extrarinse = sum(horzcat(rinse.extraall.(n{i}).s.*fsgm.*water.washer.(n{i}).*freq.cap.(n{i}), temp_extrarinse),2); % Per capita daily extra rinse water consumption
					end
					temp_extrarinse_kg = [];
					for i = 1:length(n)
						temp_extrarinse_kg = sum(horzcat((rinse.extraall.(n{i}).s.*fsgm.*water.washer.(n{i})./mass.washer.(n{i})(:,1)).*frac.(n{i}), temp_extrarinse_kg),2); % Specific extra rinse water consumption
					end
			
					temp_extrarinse_ld = [];
					for i = 1:length(n)
						temp_extrarinse_ld = sum(horzcat(rinse.extraall.(n{i}).s.*fsgm.*water.washer.(n{i}).*frac.(n{i}), temp_extrarinse_ld),2);  % Per load extra rinse water consumption 
					end	
				else	
					temp_extrarinse = [];
					for i = 1:length(n)
						temp_extrarinse = sum(horzcat(extrarinse.washer.(n{i}).*fsgm.*water.washer.(n{i}).*freq.cap.(n{i}), temp_extrarinse),2); % Per capita daily extra rinse water consumption
					end
					temp_extrarinse_kg = [];
					for i = 1:length(n)
						temp_extrarinse_kg = sum(horzcat((extrarinse.washer.(n{i}).*fsgm.*water.washer.(n{i})./mass.washer.(n{i})(:,1)).*frac.(n{i}), temp_extrarinse_kg),2); % Specific extra rinse water consumption
					end	
						temp_extrarinse_ld = [];
					for i = 1:length(n)
						temp_extrarinse_ld = sum(horzcat(extrarinse.washer.(n{i}).*fsgm.*water.washer.(n{i}).*frac.(n{i}), temp_extrarinse_ld),2); % Per load extra rinse water consumption
					end	
				end
			else
				temp_extrarinse = [];
				for i = 1:length(n)
					temp_extrarinse = sum(horzcat(water.lab.tot.(n{i}).*rinse.extra.(n{i}).s.*freq.cap.(n{i}), temp_extrarinse),2); % Per capita daily extra rinse water consumption for case study building
				end
				temp_extrarinse_kg = [];
				for i = 1:length(n)
					temp_extrarinse_kg = sum(horzcat((water.lab.tot.(n{i}).*rinse.extra.(n{i}).s./mass.washer.(n{i})(:,1)).*frac.(n{i}), temp_extrarinse_kg),2); % Specific extra rinse water consumption for case study building
				end
				temp_extrarinse_ld = [];
				for i = 1:length(n)
					temp_extrarinse_ld = sum(horzcat(water.lab.tot.(n{i}).*rinse.extra.(n{i}).s.*frac.(n{i}), temp_extrarinse_ld),2); % Per load extra rinse water consumption for case study building
				end	
			end
			
			%% Washer hot water energy
			for i = 1:length(n) 
				for s = 1:length(x)
					hotwaterbase.energy.washer.(n{i}).(x{s}) = hotwater.washer.(n{i}).(x{s}).*4.186.*(1/3600).*(storage_temp.*s_46 - mains_temp.*s_50); % Hot water energy consumption
				end
			end
			
			if scenario == 14
				for i = 1:length(n)
					temp_hotwashenergybase.(n{i}) = [];
					for s = 1:length(x)
						temp_hotwashenergybase.(n{i}) = sum(horzcat(hotwaterbase.energy.washer.(n{i}).(x{s}), temp_hotwashenergybase.(n{i})),2).*(1-spraytechwater); % Scenario #14 - Front-loading washer with spray technology
					end
				end
			else
				for i = 1:length(n)
					temp_hotwashenergybase.(n{i}) = [];
					for s = 1:length(x)
						temp_hotwashenergybase.(n{i}) = sum(horzcat(hotwaterbase.energy.washer.(n{i}).(x{s}), temp_hotwashenergybase.(n{i})),2); % Per load washer hot water energy consumption
					end
				end	
			end
			
			temp_hotwashenergyfbase = [];
			for i = 1:length(n)
				temp_hotwashenergyfbase = sum(horzcat(temp_hotwashenergybase.(n{i}).*freq.cap.(n{i}),temp_hotwashenergyfbase),2); % Per capita daily washer hot water energy consumption
			end
				temp_hotwashenergyfbase_kg = [];
			for i = 1:length(n)
				temp_hotwashenergyfbase_kg = sum(horzcat((temp_hotwashenergybase.(n{i})./mass.washer.(n{i})(:,1)).*frac.(n{i}),temp_hotwashenergyfbase_kg),2); % Specific washer hot water energy consumption
			end
			
			%% Building hot water system
			for w = 1:length(e)	
				for i = 1:length(n)
					for s = 1:length(x)
            			sntset = 58:58; % Building hot water system efficiency sensitivity analysis
            			absval = std(dhw.(e{w}));
            			[ndiff, setdi, ingroup] = prepdhw (w, e, sntset);
            			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
						if senseindex == 0
							dhwh.(e{w})(:,1) = c2*dhw.(e{w}) + d2;
						else
							if c2 > 1
								dhwh.(e{w})(:,1) = max(dhw.(e{w}));
							elseif c2 < 1
								dhwh.(e{w})(:,1) = min(dhw.(e{w}));
							else c2 = 1;
								dhwh.(e{w})(:,1) = dhw.(e{w});
							end
						end
					end
				end
			end
			
			for w = 1:length(e)
				for i = 1:length(n)
					for s = 1:length(x)
						hotwaterbase.bulk.washer.(e{w}).(n{i}).(x{s}) = hotwater.washer.(n{i}).(x{s}).*4.186.*(1/3600).*(bldg_temp.*s_47 - mains_temp.*s_50).*dhwh.(e{w}); % Building hot water system energy consumption
					end
				end
			end
			
			for w = 1:length(e)
				for i = 1:length(n)
					temp_dhw.bulkbase.(e{w}).(n{i}) = [];
					for s = 1:length(x)
						temp_dhw.bulkbase.(e{w}).(n{i}) = sum(horzcat(hotwaterbase.bulk.washer.(e{w}).(n{i}).(x{s}), temp_dhw.bulkbase.(e{w}).(n{i})),2);
					end
				end
					temp_dhwf.bulkbase.(e{w}) = [];
					for i = 1:length(n)
						temp_dhwf.bulkbase.(e{w}) = sum(horzcat(temp_dhw.bulkbase.(e{w}).(n{i}).*freq.cap.(n{i}),temp_dhwf.bulkbase.(e{w})),2); % Per capita daily building hot water system energy consumption
					end
			end
			
			%% Washer operational electricity consumption
			for i = 1:length(n)
				for kg = 1:length(lb)
        			loadbinpreelect.(n{i}).(lb{kg}) = ((mass.washer.(n{i}) >= kg) & (mass.washer.(n{i}) < (kg+1))).*kg.*electbin.(n{i}).(lb{kg});
				end
				temp_loadbinpreelect.(n{i}) = [];
				for kg = 1:length(lb)
        			temp_loadbinpreelect.(n{i}) = sum(horzcat(loadbinpreelect.(n{i}).(lb{kg}), temp_loadbinpreelect.(n{i})),2);
				end
    			for m5 = 1:length(zr1) % Washer operational electricity consumption sensitivity analysis    
        			sntset = 45:45;
        			absval = std(temp_loadbinpreelect.(n{i}));
        			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
        			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
        			if senseindex == 0
            			euro.washer.(n{i})(:,1) = c2*temp_loadbinpreelect.(n{i}) + d2;
					else
						if c2 > 1
							euro.washer.(n{i})(:,1) = max(temp_loadbinpreelect.(n{i}));
						elseif c2 < 1
							euro.washer.(n{i})(:,1) = min(temp_loadbinpreelect.(n{i}));
						else c2 = 1;
							euro.washer.(n{i})(:,1) = temp_loadbinpreelect.(n{i});
						end
        			end
    			end
			end
			
				temp_opereurof = [];
				for i = 1:length(n)
					temp_opereurof = sum(horzcat(euro.washer.(n{i}).*freq.cap.(n{i}),temp_opereurof),2); % Per capita daily washer operational electricity consumption
				end
			
				temp_opereurof_kg = [];
				for i = 1:length(n)
					temp_opereurof_kg = sum(horzcat((euro.washer.(n{i})./mass.washer.(n{i})(:,1)).*frac.(n{i}),temp_opereurof_kg),2); % Specific washer operational electricity consumption
				end
				
			%% Washer supplemental heating electricity consumption
			for i = 1:2 
				for s = 1:2 
					sntset = 40:41; % Washer supplemental heating energy consumption sensitivity analysis
					absval = std(hotwater.washer.(n{i}).(x{s}).*4.186.*(1/3600).*(tubtempsupp.(x{s}).s - tubtemp.(x{3}).s));
					[ndiff, setdi, ingroup] = preptemp (s, x, sntset);
					[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
					if senseindex == 0
						hotwaterbasesupp.energy.washer.(n{i}).(x{s}) = c2.*(hotwater.washer.(n{i}).(x{s}).*4.186.*(1/3600).*(tubtempsupp.(x{s}).s - tubtemp.(x{3}).s)) + d2;
					else
						if c2 > 1
							hotwaterbasesupp.energy.washer.(n{i}).(x{s})(:,1) = max(hotwater.washer.(n{i}).(x{s}).*4.186.*(1/3600).*(tubtempsupp.(x{s}).s - tubtemp.(x{3}).s));
						elseif c2 < 1
							hotwaterbasesupp.energy.washer.(n{i}).(x{s})(:,1) = min(hotwater.washer.(n{i}).(x{s}).*4.186.*(1/3600).*(tubtempsupp.(x{s}).s - tubtemp.(x{3}).s));
						else c2 = 1;
							hotwaterbasesupp.energy.washer.(n{i}).(x{s})(:,1) = hotwater.washer.(n{i}).(x{s}).*4.186.*(1/3600).*(tubtempsupp.(x{s}).s - tubtemp.(x{3}).s);
						end
					end
				end
			end
			
			for i = 1:2 
				supp_euro.(n{i}) = [];
				for s = 1:2
					supp_euro.(n{i}) = sum(horzcat(hotwaterbasesupp.energy.washer.(n{i}).(x{s}), supp_euro.(n{i})),2); % Per load washer supplemental heating energy consumption
				end
			end
				
			supp_eurof = [];
			for i = 1:2
				supp_eurof = sum(horzcat(supp_euro.(n{i}).*freq.cap.(n{i}),supp_eurof),2); % Per capita daily washer supplemental heating energy consumption
			end
			
			supp_eurof_kg = [];
			for i = 1:2
				supp_eurof_kg = sum(horzcat((supp_euro.(n{i})./mass.washer.(n{i})(:,1)).*frac.(n{i}),supp_eurof_kg),2); % Specific washer supplemental heating energy consumption
			end
			
			suppheating = supp_eurof;
			
			%% Dryer electricity consumption
			for i = 1:length(n)
    			sntset = 31:34; % Washer final spin speed sensitivity analysis
    			absval = std(mainrpm.(n{i}));
    			[ndiff, setdi, ingroup] = preppgrm (i, n, sntset);
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
				if senseindex == 0
        			mainrpmsensi.(n{i})(:,1) = c2*mainrpm.(n{i}) + d2;
				else
					if c2 > 1
						mainrpmsensi.(n{i})(:,1) = max(mainrpm.(n{i}));
					elseif c2 < 1
						mainrpmsensi.(n{i})(:,1) = min(mainrpm.(n{i}));
					else c2 = 1;
						mainrpmsensi.(n{i})(:,1) = mainrpm.(n{i});
					end
				end
	
				for u = 1:length(y)
					for i2 = 1:length(n2) % RMC of the wet load
						moistbinpre.(n{i}).(n2{i2}) = ((mainrpmsensi.(n{i})(:,1) >= n2r(i2)) & (mainrpmsensi.(n{i})(:,1) < (n2r(i2)+200))).*moistbin.(n{i}).(n2{i2});
					end	
					
					temp_moistbinpre.(n{i}) = [];
					for i2 = 1:length(n2)
						temp_moistbinpre.(n{i}) = sum(horzcat(moistbinpre.(n{i}).(n2{i2}), temp_moistbinpre.(n{i})),2);
					end		
					if whichmodel == 1
						rem.washer.(n{i}) = temp_moistbinpre.(n{i}); % RMC of the wet load
					else
						rem.washer.(n{i}) = rmcwetlab.(n{i}); % RMC of the wet load for the case study building
					end
					rmc.washer.(n{i}) = mass.washer.(n{i}).*rem.washer.(n{i}); % RMC (mass) remaining in the clothes after washer final spin
   			
        			sntset = 67:69; % RMC of the dry load sensitivity analysis
        			absval = std(rema.dryer.lit.(y{u}));
        			[ndiff, setdi, ingroup] = prepdry (u, y, sntset);
        			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
					if senseindex == 0
						rem.dryer.(n{i}).(y{u})(:,1) = c2*rema.dryer.lit.(y{u}) + d2;
					else 
						if c2 > 1
							rem.dryer.(n{i}).(y{u})(:,1) = max(rema.dryer.lit.(y{u}));
						elseif c2 < 1
							rem.dryer.(n{i}).(y{u})(:,1) = min(rema.dryer.lit.(y{u}));
						else c2 = 1;
							rem.dryer.(n{i}).(y{u})(:,1) = rema.dryer.lit.(y{u});
						end
					end
					rmc.dryer.(n{i}).(y{u}) = mass.washer.(n{i}).*rem.dryer.(n{i}).(y{u}); % RMC (mass) remaining in the clothes after dried in the dryer
					rmc.net.(n{i}).(y{u}) = rmc.washer.(n{i}) - rmc.dryer.(n{i}).(y{u}); % RMC net mass of moisture 
			
					if u == 1
						for m5 = 1:length(zr1)       
                			sntset = 42:42; % Dryer automatic termination correction factor for more dry sensitivity analysis
                			absval = std(autoextrafactor);
                			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
                			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
                			if senseindex == 0
                    			autofactor(:,1) = c2*autoextrafactor + d2;
							else
								if c2 > 1
									autofactor(:,1) = max(autoextrafactor);
								elseif c2 < 1
									autofactor(:,1) = min(autoextrafactor);
								else c2 = 1;
									autofactor(:,1) = autoextrafactor;
								end
                			end
						end              
					else
						for m5 = 1:length(zr1) % Dryer automatic termination correction factor for normal dry sensitivity analysis   
                			sntset = 13:13;
                			absval = std(autonormalfactor);
                			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
                			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
							if senseindex == 0
                    			autofactor(:,1) = c2*autonormalfactor + d2; 
							else
								if c2 > 1
									autofactor(:,1) = max(autonormalfactor);
								elseif c2 < 1
									autofactor(:,1) = min(autonormalfactor);
								else c2 = 1;
									autofactor(:,1) = autonormalfactor;
								end
							end
						end
					end
				
					for i3 = 1:length(n3)
						heatprecorr.(n{i}).(y{u}).(n3{i3}) = ((mass.washer.(n{i}) >= n3mass(i3)) & (mass.washer.(n{i}) < (n3mass(i3+1)))).*n3elect(i3).*((1./n3rmc(i3)).*(rmc.net.(n{i}).(y{u})./mass.washer.(n{i}))).*rmc.net.(n{i}).(y{u}); % 'h' Dryer electricity correction factor
					end	
					
					temp_heatpre.(n{i}).(y{u}) = [];
					for i3 = 1:length(n3)
						temp_heatpre.(n{i}).(y{u}) = sum(horzcat(heatprecorr.(n{i}).(y{u}).(n3{i3}), temp_heatpre.(n{i}).(y{u})),2);
					end
						
					for m5 = 1:length(zr1)          
						sntset = 43:43; % Dryer electricity consumption sensitivity analysis
						absval = std(temp_heatpre.(n{i}).(y{u}));
						[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
						[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
						if senseindex == 0 & scenario == 12
							electricity.dryer.temp.(n{i}).(y{u})(:,1) = temp_heatpre.(n{i}).(y{u}).*(heatpumpenergysavings).*autofactor; % Scenario #20 - Heat pump clothes dryer
						elseif senseindex == 0
							electricity.dryer.temp.(n{i}).(y{u})(:,1) = (c2.*temp_heatpre.(n{i}).(y{u}) + d2).*autofactor;
						else
							if c2 > 1
								electricity.dryer.temp.(n{i}).(y{u})(:,1) = max(temp_heatpre.(n{i}).(y{u})).*autofactor;
							elseif c2 < 1
								electricity.dryer.temp.(n{i}).(y{u})(:,1) = min(temp_heatpre.(n{i}).(y{u})).*autofactor;
							else c2 = 1;
								electricity.dryer.temp.(n{i}).(y{u})(:,1) = temp_heatpre.(n{i}).(y{u}).*autofactor;
							end
						end
					end  	
				end
			end
			
			dry_main = [];
			for u = 1:length(y)
				dry_main = sum(horzcat(drymain.(y{u}), dry_main),2);
			end
			for u = 1:length(y)
				tempdrystd.(y{u}) = drymain.(y{u})./dry_main;
			end
			
			for u = 1:length(y)
				sntset = 26:28; % Fraction of loads by dryness level sensitivity analysis
				absval = std(tempdrystd.(y{u}));
				[ndiff, setdi, ingroup] = prepdry (u, y, sntset);
				[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);
				if senseindex == 0
					dryness.(y{u})(:,1) = c2*tempdrystd.(y{u}) + d2;
				else
					if c2 > 1
						dryness.(y{u})(:,1) = max(tempdrystd.(y{u}));
					elseif c2 < 1
						dryness.(y{u})(:,1) = min(tempdrystd.(y{u}));
					else c2 = 1;
						dryness.(y{u})(:,1) = tempdrystd.(y{u});
					end
				end
			end
			
			for i = 1:length(n)
				for u = 1:length(y)
        			for m5 = 1:length(zr1)        
            			sntset = 44:44; % Fraction of loads dried in the dryer sensitivity analysis
            			absval = std(dryfract);
            			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
            			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
						if senseindex == 0
							electricity.dryer.(n{i}).(y{u})(:,1) = (c2*dryfract + d2).*electricity.dryer.temp.(n{i}).(y{u}).*dryness.(y{u});
						else
							if c2 > 1
								electricity.dryer.(n{i}).(y{u})(:,1) = max(dryfract).*electricity.dryer.temp.(n{i}).(y{u}).*dryness.(y{u});
							elseif c2 < 1
								electricity.dryer.(n{i}).(y{u})(:,1) = min(dryfract).*electricity.dryer.temp.(n{i}).(y{u}).*dryness.(y{u});
							else c2 = 1;
								electricity.dryer.(n{i}).(y{u})(:,1) = dryfract.*electricity.dryer.temp.(n{i}).(y{u}).*dryness.(y{u});
							end
						end
        			end 
				end
			end
			
			for i = 1:length(n)
				temp_dryer.(n{i}) = [];
				for u = 1:length(y)
					temp_dryer.(n{i}) = sum(horzcat(electricity.dryer.(n{i}).(y{u}), temp_dryer.(n{i})),2); % Per load dryer electricity consumption
				end	
			end
			
				temp_dryerf = [];
				for i = 1:length(n)
        			temp_dryerf = sum(horzcat(temp_dryer.(n{i}).*freq.cap.(n{i}),temp_dryerf),2); % Per capita daily dryer electricity consumption
				end
			
				temp_dryerf_kg = [];	
				for i = 1:length(n)
 					temp_dryerf_kg = sum(horzcat((temp_dryer.(n{i})./mass.washer.(n{i})(:,1)).*frac.(n{i}),temp_dryerf_kg),2); % Specific dryer electricity consumption
				end
				
				temp_rpm = [];
				for i = 1:length(n)
        			temp_rpm = sum(horzcat(mainrpmsensi.(n{i}).*frac.(n{i}),temp_rpm),2); % Washer final spin speed
				end
				
				temp_rmc_wsh = [];
				for i = 1:length(n)
        			temp_rmc_wsh = sum(horzcat(rem.washer.(n{i}).*frac.(n{i}),temp_rmc_wsh),2); % RMC of wet load
				end
				
				for i = 1:length(n)
					temp_rmc_dry_pre.(n{i}) = [];
					for u = 1:length(y)
						temp_rmc_dry_pre.(n{i}) = sum(horzcat(rem.dryer.(n{i}).(y{u}).*dryness.(y{u}), temp_rmc_dry_pre.(n{i})),2);
					end
				end
				temp_rmc_dry = [];
				for i = 1:length(n)
        			temp_rmc_dry = sum(horzcat(temp_rmc_dry_pre.(n{i}).*frac.(n{i}),temp_rmc_dry),2); % RMC of dry load
				end
					
			%% Energy source-to-site factors
			for m5 = 1:length(zr1)  
    			sntset = 52:52; % Source-to-site energy factor for electricity sensitivity analysis  
    			absval = std(source.electricity.s);
    			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
        			if senseindex == 0
            			sysele(:,1) = c2*source.electricity.s + d2;
					else
						if c2 > 1
							sysele(:,1) = max(source.electricity.s);
						elseif c2 < 1
							sysele(:,1) = min(source.electricity.s);
						else c2 = 1;
							sysele(:,1) = source.electricity.s;
						end
        			end
			end   
			
			for m5 = 1:length(zr1)            
    			sntset = 51:51; % Source-to-site energy factor for natural gas sensitivity analysis
    			absval = std(source.gas.s);
    			[ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset);
    			[d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist);           
        			if senseindex == 0
            			sysgas(:,1) = c2*source.gas.s + d2;
					else
						if c2 > 1
							sysgas(:,1) = max(source.gas.s);
						elseif c2 < 1
							sysgas(:,1) = min(source.gas.s);
						else c2 = 1;
							sysgas(:,1) = source.gas.s;
						end
        			end
			end
			
			%% Building hot water recirculating pump electricity consumption
			pump = (2.*0.667.*24)./(741.*68); % Per liter pump electricity consumption. Assumptions: two pumps (rated at 0.667 kW full load amperage) operating 24 h/d, and 741 tenants in the case study building each consuming 68 L/(cap*d) hot water.
			
			%% Washing frequency
			temp_freq = [];
			for i = 1:length(n)
				temp_freq = sum(horzcat(freq.cap.(n{i}), temp_freq),2);
			end
			washingfrequency = temp_freq;
			
			dummy = s_72;
			
			%% Summary
			washerwater = temp_water + temp_extrarinse; % Washer water consumption
			
			if scenario == 14 | scenario == 21 % Scenario #14 - Front-loading washer with spray technology; Scenario #12 - Water heated internally in the washer
				washerhotwater = 0;
			else	
				washerhotwater = temp_hotwashf;
			end
			
			if scenario == 14 | scenario == 21 % Scenario #14 - Front-loading washer with spray technology; Scenario #12 - Water heated internally in the washer
				washeropereuro = temp_opereurof + suppheating + temp_hotwashenergyfbase;
			elseif scenario == 13
				washeropereuro = (temp_opereurof); % Scenario #1 - Conventional top-loading washer
			else
				washeropereuro = (temp_opereurof + suppheating); % Washer electricity consumption
			end
			
			if scenario == 14 | scenario == 21 % Scenario #14 - Front-loading washer with spray technology; Scenario #12 - Water heated internally in the washer
				washerhotenergy = 0; % Washer hot water energy consumption
			else	
				washerhotenergy = temp_hotwashenergyfbase;
			end
			
			if scenario == 14 | scenario == 21 % Scenario #14 - Front-loading washer with spray technology; Scenario #12 - Water heated internally in the washer
				washerenergy = washeropereuro;
			else	
				washerenergy = washeropereuro + washerhotenergy; % Washer energy consumption
			end
			
			dryerenergyus = temp_dryerf; % Dryer electricity consumption
			
			if scenario == 14 | scenario == 21 % Scenario #14 - Front-loading washer with spray technology; Scenario #12 - Water heated internally in the washer
				centralgascons = 0; % Building hot water system natural gas energy consumption
				centralgasloss = 0; % Building hot water system natural gas energy losses
				centralelectcons = 0; % Building hot water system pump electricity consumption
			else	
				centralgascons = temp_dhwf.bulkbase.central; % Building hot water system natural gas energy consumption
				centralgasloss = centralgascons - washerhotenergy; % Building hot water system natural gas energy losses
				centralelectcons = pump.*washerhotwater; % Building hot water system pump electricity consumption
			end
				uselectcons = washeropereuro.*sysele; % Washer electricity production
				uselectloss = uselectcons - washeropereuro; % Washer electricity production losses
			
				dryerenergycons = dryerenergyus.*sysele; % Dryer electricity production
				dryerenergyloss = dryerenergycons - dryerenergyus; % Dryer electricity production losses
				
			if scenario == 14 | scenario == 21 % Scenario #14 - Front-loading washer with spray technology; Scenario #12 - Water heated internally in the washer
				centralgassyscons = 0; % Building hot water system natural gas production
				centralgassysloss = 0; % Building hot water system natural gas losses
				centralelectsyscons = 0; % Building hot water system pump electricity consumption
				centralelectsysloss = 0; % Building water system pump electricity losses
			else	
				centralgassyscons = centralgascons.*sysgas; % Building hot water system natural gas production
				centralgassysloss = centralgassyscons - centralgascons; % Building hot water system losses
				centralelectsyscons = centralelectcons.*sysele; % Building hot water system pump electricity consumption
				centralelectsysloss = centralelectsyscons - centralelectcons; % Building hot water system pump electricity losses
			end	
				systembasewasherenergy = uselectcons + centralgassyscons; % Washer primary energy
				systembasedryerenergy = dryerenergycons; % Dryer primary energy
				systembaselaundryenergy = uselectcons + centralgassyscons + dryerenergycons; % Laundry primary energy	
				
				laundryenergy = washerenergy + dryerenergyus; % Laundry energy consumption
			
			%% Sensitivity analysis compilation
			washerwater_sensi = washerwater;
			eval(['LPD1.','T',num2str(snt),'= func(washerwater_sensi);']); % Washer water consumption
			
			washerenergy_sensi = washerenergy;
			eval(['LPD2.','T',num2str(snt),'= func(washerenergy_sensi);']); % Washer energy consumption
			
			dryerenergyus_sensi = dryerenergyus;
			eval(['LPD3.','T',num2str(snt),'= func(dryerenergyus_sensi);']); % Dryer energy consumption
			
			systembasewasherenergy_sensi = systembasewasherenergy;
			eval(['LPD4.','T',num2str(snt),'= func(systembasewasherenergy_sensi);']); % Washer primary energy consumption
			
			systembasedryerenergy_sensi = systembasedryerenergy;
			eval(['LPD5.','T',num2str(snt),'= func(systembasedryerenergy_sensi);']); % Dryer primary energy consumption
				
			%% Scenario output table
			output.scenario = scenario;
			output.washerelectricity = func(washeropereuro); % Washer electricity consumption
			output.washerelectricitylosses = func(uselectloss); % Washer electricity production losses
			output.washerhotwaterenergy = func(washerhotenergy); % Building hot water system energy consumption
			output.buildingdhwsystemlosses = func(centralgasloss); % Building hot water system energy losses
			output.gasproductionlosses = func(centralgassysloss); % System gas production and transportation losses
			output.buildingpumpelectricity = func(centralelectcons); % Building hot water system pumping electricity consumption
			output.buildingpumpelectricitylosses = func(centralelectsysloss); % Building hot water system pumping electricity production losses
			output.dryerelectricity = func(dryerenergyus); % Dryer electricity consumption
			output.dryerelectricitylosses = func(dryerenergyloss); % Dryer electricity production losses
			output.washerwater = func(washerwater); % Washer water consumption
			output.washerhotwater = func(washerhotwater); % Washer hot water consumption
			output.washingfrequency = func(temp_freq); % Washing frequency
			output.loadmass = func(temp_mass); % Load mass
			output.clothing = func(clothing_all); % Clothing
			output.washerdryer = func(washeropereuro)+func(uselectloss)+func(washerhotenergy)+func(centralgasloss)+func(centralgassysloss)+func(centralelectcons)+func(centralelectsysloss)+func(dryerenergyus)+func(dryerenergyloss); % Washer + dryer primary energy consumption
			output.rpm = func(temp_rpm); % Final spin speed
			output.rmcwsh = func(temp_rmc_wsh); % RMC wet
			output.rmcdry = func(temp_rmc_dry); % RMC dry
			output.specwshwtr = func(temp_water_kg + temp_extrarinse_kg); % Specific washer water consumption
			output.loadwshwtr = func(temp_water_ld + temp_extrarinse_ld); % Per load washer water consumption
			output.specwshnrg = func(temp_hotwashenergyfbase_kg + temp_opereurof_kg + supp_eurof_kg); % Specific washer energy consumption
			output.specdrynrg = func(temp_dryerf_kg); % Specific dryer energy consumption	
			eval(['SCN.T',num2str(scenario),' = output;']);
		end
	end
	
	%% Scenario analysis for Fig. 6. Per capita washer water and primary energy consumption...
	for scenario = [s_nameSC.Var1]'
		eval(['SCN',num2str(scenario),'_bar = structfun(func,SCN.T',num2str(scenario),');']);
	end
	
	Scenario_varnames = array2table(fieldnames(output));
	Scenario_varnames.Properties.VariableNames = {'Scenario'};
	
	Scenario_table = horzcat(SCN17_bar,	SCN6_bar,	SCN7_bar,	SCN23_bar,	SCN3_bar,	SCN25_bar,	SCN4_bar,	SCN5_bar,	SCN26_bar,	SCN30_bar,	SCN9_bar,	SCN10_bar,	SCN15_bar,	SCN13_bar,	SCN14_bar,	SCN16_bar,	SCN12_bar,	SCN21_bar,	SCN22_bar, SCN31_bar, SCN0_bar);
	Scenario_presort = array2table(Scenario_table');
	Scenario_presortjoin = horzcat(s_nameSC.scenario,Scenario_presort);
	Scenario_presortjoin2 = sortrows(Scenario_presortjoin,17,'ascend');
	Scenario_headers = Scenario_presortjoin2.Var1_1';
	Scenario_headersinsert={'scenario_name',Scenario_headers{1:end}};
	Scenario_sort = array2table(sortrows(Scenario_table',16,'ascend')');
	Scenario_output = horzcat(Scenario_varnames,Scenario_sort);
	Scenario_output2 = table2cell(Scenario_output);
	Scenario_output3 = [Scenario_headersinsert;Scenario_output2];
	
	%% Sensitivity analysis compilation
	for dlg = 1:5
		eval(['LPD',num2str(dlg),'_bar = structfun(func,LPD',num2str(dlg),');']);
	end
	
	%% Metered and modeled data for Fig. 5. Cumulative distribution functions...
	metered.(g).water = metered.(g).(outlier).water; % Washer water consumption
	actual = prep_actual(metered.(g).water);
	modif = prep_mod(washerwater);
	PBIAS.total = pbias(actual, modif);
	NSE.total = nse(actual, modif);
	subplot (2,2,1);
	h(1,1) = cdfplot(actual);
	hold on;
	h(1,2) = cdfplot(modif);
	set(gca, 'FontName', 'Times New Roman');
	title({'Per capita washer','water consumption'});
	xlabel('L/(cap � d)')
	ylabel('Cumulative probability')
	yticks([0:0.25:1]);
	set( h(:,1), 'LineStyle', '--', 'Color', 'b');
	
	if whichmodel == 1
		set( h(:,2), 'LineStyle', '-', 'Color', 'k');
		legend({'Metered in case study building','Modeled based on literature input data'},'Location','best','Orientation','vertical'); % Legend
	else
		set( h(:,2), 'LineStyle', '-', 'Color', 'b');
		legend({'','Modeled based on literature input data','Metered in case study building','Modeled based on case study building and literature input data'},'Location','best','Orientation','vertical')
	
		labels = {'a)', 'b)', 'c)', 'd)'}; % Fig5 panel labels
		for k = 1:4
    		subplot(2,2,k)  	
    		xLimits = xlim; % Get axis limits
    		yLimits = ylim;   	
    		xText = xLimits(2) - 0.05 * range(xLimits); % Position text near top-right with a small offset
    		yText = yLimits(2) - 0.05 * range(yLimits);
		
    		text(xText, yText, labels{k}, ...
         		'FontWeight', 'bold', 'FontSize', 12, ...
		 		'FontName', 'Times New Roman', ...
         		'HorizontalAlignment', 'right', ...
         		'VerticalAlignment', 'top')
		end
	
	fig = gcf;  % get current figure handle
	fig.Name = 'Fig5'; % rename figure
	fig.NumberTitle = 'off';
	
	end
	xticks([0:10:80]);
	
	metered.(g).hotwater = metered.(g).(outlier).hotwater; % Hot water consumption
	actual = prep_actual(metered.(g).hotwater);
	modif = prep_mod(washerhotwater);
	PBIAS.hot = pbias(actual, modif);
	NSE.hot = nse(actual, modif);
	subplot (2,2,2);
	h(1,1) = cdfplot(actual);
	hold on;
	h(1,2) = cdfplot(modif);
	set(gca, 'FontName', 'Times New Roman');
	title({'Per capita building','hot water consumption'});
	xlabel('L/(cap � d)')
	ylabel('Cumulative probability')
	yticks([0:0.25:1]);
	set( h(:,1), 'LineStyle', '--', 'Color', 'b');
	if whichmodel == 1
		set( h(:,2), 'LineStyle', '-', 'Color', 'k');
	else
		set( h(:,2), 'LineStyle', '-', 'Color', 'b');
	end
	xticks([0:1:7]);
	
	metered.(g).washerenergy = metered.(g).(outlier).washerenergy; % Washer energy consumption
	actual = prep_actual(metered.(g).washerenergy);
	modif = prep_mod(washerenergy);
	PBIAS.washerelect = pbias(actual, modif);
	NSE.washerelect = nse(actual, modif);
	if whichmodel==1
		subplot (2,2,3);
		h(1,1) = cdfplot(actual);
		hold on;
		h(1,2) = cdfplot(modif);
		set(gca, 'FontName', 'Times New Roman');
		title({'Per capita washer','energy consumption'});
		xlabel('kWh/(cap � d)')
		ylabel('Cumulative probability')
		yticks([0:0.25:1]);
		set( h(:,1), 'LineStyle', '--', 'Color', 'b');
		set( h(:,2), 'LineStyle', '-', 'Color', 'k');
		xticks([0:0.1:0.6]);
	end
	
	metered.(g).dryerelectricity = metered.(g).(outlier).dryerelectricity; % Dryer electricity consumption
	actual = prep_actual(metered.(g).dryerelectricity);
	modif = prep_mod(dryerenergyus);
	PBIAS.electdryer = pbias(actual, modif);
	NSE.electdryer = nse(actual, modif);
		subplot (2,2,4);
		h(1,1) = cdfplot(actual);
		hold on;
		h(1,2) = cdfplot(modif);
		set(gca, 'FontName', 'Times New Roman');
		title({'Per capita dryer','electricity consumption'});
		xlabel('kWh/(cap � d)')
		ylabel('Cumulative probability')
		yticks([0:0.25:1]);
		set( h(:,1), 'LineStyle', '--', 'Color', 'b');
	if whichmodel == 1
		set( h(:,2), 'LineStyle', '-', 'Color', 'k');
	else
		set( h(:,2), 'LineStyle', '-', 'Color', 'b');
	end
	xticks([0:0.5:4]);
	
	stat = {'nanmedian' 'nanmean' 'min' 'max'};
	stream2 = {'washerwater' 'washerhotwater'  'washerenergy' 'dryerenergyus'}; % Based on literature model
	stream = {'water' 'hotwater' 'washerenergy' 'dryerelectricity'}; % Based on case study building data + literature data model
	for st3 = 1:length(stat)
		for st4 = 1:length(stream)
			for st5 = 1:length(stream2)
				eval(['metrics.metered.',(stat{st3}),'.',(stream{st4}),'=',(stat{st3}),'(metered.(g).(outlier).',(stream{st4}),');']);
				eval(['metrics.modeled.',(stat{st3}),'.',(stream2{st5}),'=',(stat{st3}),'(',(stream2{st5}),');']);
			end	
		end
	end
	
	% Table 1. Metered and modeled washer water and hot water consumption...
	sd =3;
	sg = 'decimals';
	Table1 = table(round(structfun(func,metrics.metered.nanmean),sd,sg), round(structfun(func,metrics.modeled.nanmedian),sd,sg), round(structfun(func,metrics.modeled.min),sd,sg), round(structfun(func,metrics.modeled.max),sd,sg), round(structfun(func,PBIAS),sd,sg), round(structfun(func,NSE),sd,sg), 'RowNames', stream);
	Table1.Properties.VariableNames = {'Metered_weighted_average','Modeled_median','Modeled_minimum', 'Modeled_maximum','PBIAS','NSE'};
	
	if whichmodel == 0
		Table1.Modeled_median(3:3)="#N/A";
		Table1.Modeled_minimum(3:3)="#N/A";
		Table1.Modeled_maximum(3:3)="#N/A";
		Table1.PBIAS(3:3)="#N/A";
		Table1.NSE(3:3)="#N/A";
	end
	
	%  SM E.1. Sensitivity index for the model input variables.
	Table2 = table(s_name.Var1, s_name.central, LPD1_bar(:,1), LPD2_bar(:,1), LPD3_bar(:,1), LPD4_bar(:,1), LPD5_bar(:,1));
	Table2.Properties.VariableNames = {'SN', 'Input', 'Washer_water', 'Washer_energy', 'Dryer_energy','System_washer_energy' 'System_dryer_energy'};
	
			if whichmodel == 1 && whichsensi == 0
				Table1_lit_model_performance = Table1; % Table 1. Metered and nodeled (literature) washer water and hot water consumption, washer and dryer energy consumption, and performance metrics.
				openvar Table1_lit_model_performance
				
				Fig6_scenario_table = Scenario_output3; % Scenario analysis data for Fig. 6. Per capita washer water and primary energy consumption...
				openvar Fig6_scenario_table
			
			elseif whichmodel == 0 && whichsensi == 0
				Table1_lit_bldg_model_performance = Table1; % Table 1. Metered and modeled (literature + case study building) washer water and hot water consumption, washer and dryer energy consumption, and performance metrics.
				openvar Table1_lit_bldg_model_performance
			
			elseif whichmodel == 1 && whichsensi == 1 % Max for determination of SI index
				SImax_table = Table2;
			
			elseif whichmodel == 1 && whichsensi == 2 % Min for determination of SI index
				SImin_table = Table2;
			else
			
			end
		end

	% SI index table for SM E.1. Sensitivity index...
	SImin_values = SImin_table{:, 3:7}; % Extract numeric data from columns 3 to 7
	SImax_values = SImax_table{:, 3:7};
	min_vals = min(SImin_values, SImax_values); % Element-wise min and max
	max_vals = max(SImax_values, SImin_values);
	max_vals(max_vals == 0) = NaN; % Avoid division by zero
	SIindex = 1 - (min_vals ./ max_vals); % Calculate sensitivity index
	SIindex_table = array2table(SIindex, ...
    'VariableNames', SImin_table.Properties.VariableNames(3:7)); % Convert result to table with original variable names
	SIindex_table = [SImin_table(:, 1:2), SIindex_table]; 	% Add text columns (e.g., names) back in
	openvar SIindex_table
end

clc
fprintf('%s completed successfully',name)

beep

%% Functions
function b = prep_actual(v) % Preparing actual data
actual = sort(v);
actual(isnan(actual))=[];
b = actual;
end

function c = prep_mod(v) % Preparing modeled data
mod = sort(v);
mod(isnan(mod))=[];
c = mod;
end

function d = pbias(actual, modif) % Preparing and analyzing PBIAS performance metrics data
actual_rank = (tiedrank(actual) / length(actual))*100;
mod_rank = prctile(modif,actual_rank);
d = cell2mat({(sum(actual - mod_rank).*100)/sum(actual)});
end

function e = nse(actual, modif) % Preparing and analyzing NSE performance metrics data
actual_rank = (tiedrank(actual) / length(actual))*100;
mod_rank = prctile(modif,actual_rank);
e = cell2mat({1 - (sum((actual - mod_rank).^2)/sum((actual - mean(actual)).^2))});
end

function [ndiff, setdi, ingroup] = preppgrm (i, n, sntset) % Preparing program data
ndiff = 1:length(n);
setdi = intersect(i,ndiff);
ingroup = sntset(i);
end

function [ndiff, setdi, ingroup] = prephhs (htx, hx, sntset)  % Preparing household size data
ndiff = 1:length(hx);
setdi = intersect(htx,ndiff);
ingroup = sntset(htx);
end

function [ndiff, setdi, ingroup] = prepsensi (m5, zr1, sntset)  % Preparing sensitivity analysis data
ndiff = 1:length(zr1);
setdi = intersect(m5,ndiff);
ingroup = sntset(m5);
end

function [ndiff, setdi, ingroup] = preptemp (s, x, sntset) % Preparing wash water temperature data
ndiff = 1:length(x);
setdi = intersect(s,ndiff);
ingroup = sntset(s);
end

function [ndiff, setdi, ingroup] = prepdry (u, y, sntset)  % Preparing dryer data
ndiff = 1:length(y);
setdi = intersect(u,ndiff);
ingroup = sntset(u);
end

function [ndiff, setdi, ingroup] = prepdhw (w, e, sntset) % Preparing hot water system data
ndiff = 1:length(e);
setdi = intersect(w,ndiff);
ingroup = sntset(w);
end

function [d2, c2] = sensifunc (ndiff, setdi, ingroup, snt, sntset, sensibase, sensichange, absval, sensitype, redist) % Sensitivity analysis
setdnd = setdiff(ndiff,setdi);
sensiweight = sensibase - ((sensichange - sensibase)/length(setdnd));	
absvalweight = absval/length(setdnd);
	if ismember(snt,sntset)
		if snt == ingroup
			if sensitype == 0											
			d2 = absval;													
			c2 = sensibase;													
			else														
			d2 = 0;															
			c2 = sensichange;													
			end													
		else
			if redist == 0																				
				d2 = 0;												
				c2 = sensibase;
			else															
				if sensitype == 0
				d2 = -absvalweight;
				c2 = sensibase;
				else
				d2 = 0;
				c2 = sensiweight;	

				end
			end			
		end
	else
	d2 = 0;
	c2 = sensibase;
	end
end

function [distr] = extdist(hh1, hh2, hhs, ss)  % Preparing washing frequency of use data
[pdca,gn,gl] = fitdist(hh1,'Normal','By',hh2);
dist1 = pdca{1, hhs}.InputData.data; 
dist2 = random(truncate(fitdist(dist1,'Normal'),min(dist1),max(dist1)),ss,1);
distr = datasample(dist2,ss,1);
end